Whittle–Matérn fields with general smoothness

Introduction

In this vignette we will introduce how to fit Whittle–Matérn fields with general smoothness based on finite element and rational approximations. The theory for this approach is provided in Bolin et al. (2022) and Xiong, Simas, and Bolin (2022). For the implementation, make use of the rSPDE package for the rational approximations.

These models are thus implemented using finite element approximations. Such approximations are not needed in for integer smoothness parameters, and for the details about the exact models we refer to the vignettes

• Whittle–Matérn fields on metric graphs
• INLA and inlabru interfaces

For further details on the construction of metric graphs, see Working with metric graphs

Constructing the graph and the mesh

We begin by loading the rSPDE and MetricGraph packages:

  library(rSPDE)
  library(MetricGraph)

As an example, we consider the following metric graph

  library(sp)
  line1 <- Line(rbind(c(0,0),c(1,0)))
  line2 <- Line(rbind(c(0,0),c(0,1)))
  line3 <- Line(rbind(c(0,1),c(-1,1)))
  theta <- seq(from=pi,to=3*pi/2,length.out = 20)
  line4 <- Line(cbind(sin(theta),1+ cos(theta)))
  Lines = SpatialLines(list(Lines(list(line1),ID="1"),
                                Lines(list(line2),ID="2"),
                                Lines(list(line3),ID="3"),
                                Lines(list(line4),ID="4")))
  graph <- metric_graph$new(lines = Lines)
  graph$plot()

To construct a FEM approximation of a Whittle–Matérn field with general smoothness, we must first construct a mesh on the graph.

  graph$build_mesh(h = 0.5)
  graph$plot(mesh=TRUE)

In the command build_mesh, the argument h decides the largest spacing between nodes in the mesh. As can be seen in the plot, the mesh is very coarse, so let’s reduce the value of h and rebuild the mesh:

graph$build_mesh(h = 0.01)

The next step is to build the mass and stiffness matrices for the FEM basis.

graph$compute_fem()

We are now ready to specify the model \[ (\kappa^2 - \Delta)^{\alpha} \tau u = \mathcal{W} \] for the Whittle–Matérn field \(u\). For this, we use the matern.operators function from the rSPDE package:

  sigma <- 1.3
  r <- 0.2 # range
  nu <- 0.9
  rspde.order <- 2
  kappa <- sqrt(8 * nu) / r
  op <- matern.operators(C= graph$mesh$C, G = graph$mesh$G, d = 1, 
                         nu = nu, kappa = kappa, sigma = sigma, 
                         m = rspde.order)

As can be seen in the code, we specify \(\kappa\) via the practical correlation range \(\sqrt{8\nu}/\kappa\). Also, the model is not parametrized by \(\tau, \alpha\) but instead by \(\sigma, \nu\). Here, sigma denotes the standard deviation of the field and nu is the smoothness parameter, which is related to \(\alpha\) via the relation \(\alpha = \nu + 1/2\). The object op_cov contains the matrices needed for evaluating the distribution of the stochastic weights in the FEM approximation.

Let us simulate the field \(u\) at the mesh locations and plot the result:

u <- simulate(op)
graph$plot_function(u, plotly = TRUE)

If we want to evaluate \(u(s)\) at some locations \(s_1,\ldots, s_n\), we need to multiply the weights with the FEM basis functions \(\varphi_i(s)\) evaluated at the locations. For this, we can construct the observation matrix \(\boldsymbol{\mathrm{A}}\), with elements \(A_{ij} = \varphi_j(s_i)\), which links the FEM basis functions to the locations. This can be done by the function mesh_A in the metric graph object. To illustrate this, let us simulate some observation locations on the graph and construct the matrix:

obs.per.edge <- 50
obs.loc <- NULL
for(i in 1:graph$nE) {
  obs.loc <- rbind(obs.loc,
                   cbind(rep(i,obs.per.edge), runif(obs.per.edge)))
}
n.obs <- obs.per.edge*graph$nE
A <- graph$mesh_A(obs.loc)

In the code, we generate \(50\) observation locations per edge in the graph, drawn at random. It can be noted that we assume that the observation locations are given in the format \((e, d)\) where \(e\) denotes the edge of the observation and \(d\) is the position on the edge, i.e., the relative distance from the first vertex of the edge.

To compute the precision matrix from the covariance-based rational approximation one can use the precision() method on object returned by the matern.operators() function:

  Q <- precision(op)

As an illustration of the model, let us compute the covariance function between the process at the mid point of the second edge and all other locations in the mesh. The covariances can be calculated as \[ \overline{\boldsymbol{\mathrm{A}}} \boldsymbol{\mathrm{Q}}^{-1}\overline{\boldsymbol{\mathrm{v}}}. \] Here, \(\boldsymbol{\mathrm{Q}}\) is the precision matrix obtained from the rational approximation, \(\boldsymbol{\mathrm{A}}\) is an identity matrix since we are evaluating the approximation in the nodes of the FEM mesh, \(\overline{\boldsymbol{\mathrm{v}}}\) is the \((m+1)\)-fold vertical concatenation of the vector \(\boldsymbol{\mathrm{v}}\), where \(\boldsymbol{\mathrm{v}}\) is a vector with all basis functions evaluated in \(s=0.5\).

  v <- t(graph$mesh_A(matrix(c(2,0.5),1,2)))
  v_bar <- kronecker(matrix(1, nrow = rspde.order + 1),
                     v)
  A_bar <- kronecker(matrix(1, ncol = rspde.order + 1), 
                     Matrix::Diagonal(dim(graph$mesh$V)[1]))
  c_cov <- (A_bar) %*% solve(Q, v_bar)
  graph$plot_function(c_cov, plotly = TRUE)

Using the model for inference

There is built-in support for computing log-likelihood functions and performing kriging prediction in the rSPDE package which we can use for the graph model. To illustrate this, we use the simulation to create some noisy observations of the process. We generate the observations as \(Y_i = u(s_i) + \varepsilon_i\), where \(\varepsilon_i \sim N(0,\sigma_e^2)\) is Gaussian measurement noise.

    sigma.e <- 0.1
    Y <- as.vector(A %*% u + sigma.e * rnorm(n.obs))

Let us now fit the model. To this end we first must compute the loglikelihood function as function of the parameters we want to estimate. We define the loglikelihood function parametrized using the logarithm of each parameter to avoid constrained optimization.

    mlik <- rSPDE.construct.matern.loglike(op, Y=Y, A=A, parameterization = "matern")

We now set some suitable initial values for the optimization and fit the model using optim.

  theta0 <- graph_starting_values(graph, range_par = TRUE,
                                    nu = TRUE, like_format = TRUE,
                                    manual_data = Y, log = TRUE)

  theta <- optim(theta0, mlik, method = "L-BFGS-B")
    
  print(data.frame(sigma = c(sigma, exp(theta$par[1])), 
                   range = c(r, exp(theta$par[2])),
                   nu = c(nu, exp(theta$par[3])),
                   sigma.e = c(sigma.e, exp(theta$par[4])),
                   row.names = c("Truth", "Estimates")))

##              sigma    range        nu   sigma.e
## Truth     1.300000 0.200000 0.9000000 0.1000000
## Estimates 1.385231 0.243925 0.8771482 0.1047377

  sigma_est <- exp(theta$par[1])
  kappa_est <- exp(theta$par[2])
  nu_est <- exp(theta$par[3])

  op_cov <- update(op,
                   user_kappa = kappa_est,
                   user_sigma = sigma_est,
                   user_nu = nu_est)

  u.krig <- predict(op, A = A, Aprd = Matrix::Diagonal(op$sizeC), Y = Y, sigma.e = sigma.e)

  graph$plot_function(as.vector(u.krig[[1]]), plotly = TRUE)  

Fitting a model with replicates

Let us now illustrate how to simulate a data set with replicates and then fit a model to such data. To simulate a latent model with replicates, all we do is set the nsim argument to the number of replicates.

  n.rep <- 30
  u.rep <- simulate(op, nsim = n.rep)

Now, let us generate the observed values \(Y\):

  sigma.e <- 0.3
  Y.rep <- A %*% u.rep + sigma.e * matrix(rnorm(n.obs * n.rep), ncol = n.rep)

Note that \(Y\) is a matrix with 20 columns, each column containing one replicate. Now, the remaining of the code is identical to the previous case. The rSPDE.matern.loglike() function automatically identifies the replicates from the fact that \(Y\) is a matrix with more than one column.

  theta0 <- graph_starting_values(graph, range_par = TRUE,
                                    nu = TRUE, like_format = TRUE,
                                    manual_data = Y.rep, log = TRUE)

    mlik <- rSPDE.construct.matern.loglike(op_cov, Y=Y.rep, A=A, parameterization = "matern")

    theta <- optim(theta0, mlik, method = "L-BFGS-B")
    
    print(data.frame(sigma = c(sigma, exp(theta$par[1])),
                     range = c(r, exp(theta$par[2])),
                     nu = c(nu, exp(theta$par[3])),
                     sigma.e = c(sigma.e, exp(theta$par[4])),
                     row.names = c("Truth", "Estimates")))

##              sigma     range        nu   sigma.e
## Truth     1.300000 0.2000000 0.9000000 0.3000000
## Estimates 1.275934 0.1898274 0.8898979 0.3052531

Using the R-INLA implementation

We also have an R-INLA implementation of the rational SPDE approach for metric graphs.

We begin by defining the model by using the rspde.metric_graph() function. This function contains the same arguments as the function rspde.matern(). We refer the reader to the R-INLA implementation of the rational SPDE approach vignette for further details.

We begin by adding the observations to the graph:

graph$add_observations(data = data.frame(y=Y, 
                          edge_number = obs.loc[,1], 
                          distance_on_edge = obs.loc[,2]), 
                          normalized = TRUE)

Let us create the model object:

  library(INLA)
  rspde_model <- rspde.metric_graph(graph)

By default, the order of the rational approximation is 2.

We can now create the A matrix and the index object by using the same functions as for the rational SPDE approach, namely, rspde.make.A() and rspde.make.index(), supplying the graph object as the mesh argument. Here we create the A matrix:

  Abar <- rspde.make.A(graph)

Now, let us create the index object:

  rspde.index <- rspde.make.index(name="field", mesh=graph)

The remaining is standard: we create the formula object, the stack object, and then fit the model by using the inla() function. So, first we create the formula object:

  f.s <- y ~ -1 + Intercept + f(field, model = rspde_model)

Now we create the inla.stack object:

  stk.dat <- inla.stack(
    data = graph_data_spde(rspde_model), A = Abar, tag = "est",
    effects =
      c(
        rspde.index,
        list(Intercept = 1)
      )
    )

Finally, we can fit the model:

  rspde_fit <- inla(f.s, data = inla.stack.data(stk.dat),
    control.inla = list(int.strategy = "eb"),
    control.predictor = list(A = inla.stack.A(stk.dat), compute = TRUE),
              inla.mode = "experimental"
  )

We can use the same functions as the rspde fitted models in inla. For instance, we can see the results in the original scale by creating the result object:

  result_fit <- rspde.result(rspde_fit, "field", rspde_model)
  summary(result_fit)

##             mean        sd 0.025quant 0.5quant 0.975quant     mode
## std.dev 1.408070 0.1800640   1.091420 1.395020   1.797600 1.367280
## range   0.241322 0.0524461   0.153785 0.236118   0.358954 0.225903
## nu      0.902068 0.1309470   0.654682 0.898850   1.165900 0.889826

Let us compare with the true values:

  result_df <- data.frame(
    parameter = c("std.dev", "range", "nu"),
    true = c(sigma, r, nu),
    mean = c(
      result_fit$summary.std.dev$mean,
      result_fit$summary.range$mean,
      result_fit$summary.nu$mean
    ),
    mode = c(
      result_fit$summary.std.dev$mode,
      result_fit$summary.range$mode,
      result_fit$summary.nu$mode
    )
  )
  print(result_df)

##   parameter true      mean      mode
## 1   std.dev  1.3 1.4080696 1.3672847
## 2     range  0.2 0.2413216 0.2259025
## 3        nu  0.9 0.9020683 0.8898255

We can also plot the posterior marginal densities with the help of the gg_df() function:

  posterior_df_fit <- gg_df(result_fit)

  library(ggplot2)

  ggplot(posterior_df_fit) + geom_line(aes(x = x, y = y)) + 
  facet_wrap(~parameter, scales = "free") + labs(y = "Density")

  pred_loc <- graph$mesh$VtE

graph$add_observations(data = data.frame(y=rep(NA,nrow(pred_loc)), 
                                edge_number = pred_loc[,1], 
                                distance_on_edge = pred_loc[,2]), 
                                normalized = TRUE)

  Abar_prd <- rspde.make.A(graph, loc=pred_loc)

  ef.prd <- 
    c(rspde.index, list(Intercept = 1))
  stk.prd <- inla.stack(
    data = graph_data_spde(rspde_model, only_pred=TRUE),
    A = Abar_prd, tag = "prd",
    effects = ef.prd
  )
  stk.all <- inla.stack(stk.dat, stk.prd)

rspde_fitprd <- inla(f.s,
  data = inla.stack.data(stk.all),
  control.predictor = list(
    A = inla.stack.A(stk.all),
    compute = TRUE, link = 1
  ),
  control.compute = list(
    return.marginals = FALSE,
    return.marginals.predictor = FALSE
  ),
  control.inla = list(int.strategy = "eb")
)

id.prd <- inla.stack.index(stk.all, "prd")$data
m.prd <- rspde_fitprd$summary.fitted.values$mean[id.prd]

  graph$plot_function(m.prd, plotly = TRUE)  

Using R-INLA implementation to fit models with replicates

Let us begin by cloning the graph and clearing the observations on the cloned graph:

graph_rep <- graph$clone()
graph_rep$clear_observations()

We will now add the data with replicates to the graph:

graph_rep$add_observations(data = data.frame(y=as.vector(Y.rep), 
                          edge_number = rep(obs.loc[,1], n.rep), 
                          distance_on_edge = rep(obs.loc[,2], n.rep),
                          repl = rep(1:n.rep, each = n.obs)), 
                          group = "repl",
                          normalized = TRUE)

To fit the model with replicates we need to create the corresponding A matrix and index object:

Abar.rep <- rspde.make.A(
  mesh = graph_rep, index = rep(1:n.obs, times = n.rep),
  repl = rep(1:n.rep, each = n.obs)
)

mesh.index.rep <- rspde.make.index(
  name = "field", mesh = graph_rep,
  n.repl = n.rep
)

Similarly, let us create a new rspde model object:

rspde_model_rep <- rspde.metric_graph(graph_rep)

Let us now create the corresponding inla.stack object, where we set the repl argument in the function graph_data_spde to __all since we want to use all replicates:

st.dat.rep <- inla.stack(
  data = graph_data_spde(rspde_model_rep,
                            repl = "__all"),
  A = Abar.rep,
  effects = mesh.index.rep
)

Observe that we need the response variable y to be a vector. We can now create the formula object, remembering that since we gave the name argument field, when creating the index, we need to pass field.repl to the formula:

f.rep <-
  y ~ -1 + f(field,
    model = rspde_model_rep,
    replicate = field.repl
  )

We can, finally, fit the model:

rspde_fit_rep <-
  inla(f.rep,
    data = inla.stack.data(st.dat.rep),
    family = "gaussian",
    control.predictor =
      list(A = inla.stack.A(st.dat.rep)),
            inla.mode = "experimental"
  )

We can obtain the estimates in the original scale with the rspde.result() function:

  result_fit_rep <- rspde.result(rspde_fit_rep, "field", rspde_model)
  summary(result_fit_rep)

##             mean        sd 0.025quant 0.5quant 0.975quant     mode
## std.dev 1.272090 0.0270196   1.220920 1.271370    1.32702 1.269320
## range   0.187300 0.0104929   0.169241 0.186351    0.21024 0.183601
## nu      0.904973 0.0686166   0.760356 0.909706    1.02751 0.926811

Let us compare with the true values of the parameters:

  result_rep_df <- data.frame(
    parameter = c("std.dev", "range", "nu"),
    true = c(sigma, r, nu),
    mean = c(
      result_fit_rep$summary.std.dev$mean,
      result_fit_rep$summary.range$mean,
      result_fit_rep$summary.nu$mean
    ),
    mode = c(
      result_fit_rep$summary.std.dev$mode,
      result_fit_rep$summary.range$mode,
      result_fit_rep$summary.nu$mode
    )
  )
  print(result_rep_df)

##   parameter true      mean      mode
## 1   std.dev  1.3 1.2720948 1.2693174
## 2     range  0.2 0.1873000 0.1836008
## 3        nu  0.9 0.9049733 0.9268108

We can also plot the posterior marginal densities with the help of the gg_df() function:

  posterior_df_fit_rep <- gg_df(result_fit_rep)

  ggplot(posterior_df_fit_rep) + geom_line(aes(x = x, y = y)) + 
  facet_wrap(~parameter, scales = "free") + labs(y = "Density")

Using inlabru implementation

The inlabru package allows us to fit models and do kriging in a straighforward manner, without having to handle A matrices, indices nor inla.stack objects. Therefore, we suggest the reader to use this implementation when using our implementation to fit real data.

Let us clear the graph, since it contains NA observations we used for prediction, add the observations again, and create a new rSPDE model object:

graph$clear_observations()
graph$add_observations(data = data.frame(y=Y, 
                          edge_number = obs.loc[,1], 
                          distance_on_edge = obs.loc[,2]), 
                          normalized = TRUE)
rspde_model <- rspde.metric_graph(graph)

Let us now load the inlabru package and create the component (which is inlabru’s formula-like object). Since we are using the data from the graph, inlabru will also obtain the locations from the graph, thus, there is no need to provide the locations. However, we need a name for the locations for using inlabru’s predict method. Therefore, we can choose any name for the location that is not a name being used in the graph’s data. In our case we will use the name loc:

    library(inlabru)
    cmp <-
    y ~ -1 + Intercept(1) + field(loc, 
                          model = rspde_model)

Now, we can directly fit the model, where we pass the name of the location variable as the loc argument of the graph_data_spde() function:

  rspde_bru_fit <-
    bru(cmp,
        data=graph_data_spde(rspde_model, loc = "loc"),
      options=list(
      family = "gaussian",
      inla.mode = "experimental")
    )

Let us now obtain the estimates of the parameters in the original scale by using the rspde.result() function:

  result_bru_fit <- rspde.result(rspde_bru_fit, "field", rspde_model)
  summary(result_bru_fit)

##             mean        sd 0.025quant 0.5quant 0.975quant     mode
## std.dev 1.408250 0.1799740   1.091840 1.395180   1.797670 1.367240
## range   0.241471 0.0526009   0.153772 0.236222   0.359531 0.225915
## nu      0.901841 0.1315030   0.653310 0.898657   1.166680 0.889757

Let us compare with the true values of the parameters:

  result_bru_df <- data.frame(
    parameter = c("std.dev", "range", "nu"),
    true = c(sigma, r, nu),
    mean = c(
      result_bru_fit$summary.std.dev$mean,
      result_bru_fit$summary.range$mean,
      result_bru_fit$summary.nu$mean
    ),
    mode = c(
      result_bru_fit$summary.std.dev$mode,
      result_bru_fit$summary.range$mode,
      result_bru_fit$summary.nu$mode
    )
  )
  print(result_bru_df)

##   parameter true      mean      mode
## 1   std.dev  1.3 1.4082541 1.3672408
## 2     range  0.2 0.2414709 0.2259154
## 3        nu  0.9 0.9018408 0.8897565

We can also plot the posterior marginal densities with the help of the gg_df() function:

  posterior_df_bru_fit <- gg_df(result_bru_fit)

  ggplot(posterior_df_bru_fit) + geom_line(aes(x = x, y = y)) + 
  facet_wrap(~parameter, scales = "free") + labs(y = "Density")

  data_prd_list <- graph$get_mesh_locations(bru = TRUE,
                                            loc = "loc")

  field_pred <- predict(rspde_model, cmp, rspde_bru_fit, data=data_prd_list, ~field)

  plot(field_pred, plotly = TRUE)  

Using inlabru to fit models with replicates

We can also use our inlabru implementation to fit models with replicates. We will consider the same data that was generated above, where the number of replicates is 30.

For this implementation we will use the rspde_model_rep object.

We can now create the component, passing the vector with the indices of the replicates as the replicate argument. For the replicate argument we use the function graph_repl_spde() to extract the vector of replicates, where we set choose __all as the replicates, since we want all replicates.

  cmp_rep <-
    y ~ -1 + Intercept(1) + field(loc, model = rspde_model_rep,
                              replicate = graph_repl_spde(rspde_model_rep, "__all"))

Now, we are ready to fit the model. Observe that we set the repl argument to __all since we will use all replicates and the loc argument to loc, which is the name we gave to the locations.

  rspde_bru_fit_rep <-
    bru(cmp_rep,
        data=graph_data_spde(rspde_model_rep, 
                         repl = "__all",
                         loc = "loc"),
      options=list(
      family = "gaussian",
      inla.mode = "experimental")
    )

We can obtain the estimates in the original scale with the rspde.result() function:

  result_bru_fit_rep <- rspde.result(rspde_bru_fit_rep, "field", rspde_model)
  summary(result_bru_fit_rep)

##             mean        sd 0.025quant 0.5quant 0.975quant     mode
## std.dev 1.280840 0.0273472   1.227870 1.280590   1.335290 1.280150
## range   0.191240 0.0106389   0.169550 0.191665   0.211052 0.193754
## nu      0.897091 0.0672550   0.787134 0.889103   1.046570 0.860258

Let us compare with the true values of the parameters:

  result_bru_rep_df <- data.frame(
    parameter = c("std.dev", "range", "nu"),
    true = c(sigma, r, nu),
    mean = c(
      result_bru_fit_rep$summary.std.dev$mean,
      result_bru_fit_rep$summary.range$mean,
      result_bru_fit_rep$summary.nu$mean
    ),
    mode = c(
      result_bru_fit_rep$summary.std.dev$mode,
      result_bru_fit_rep$summary.range$mode,
      result_bru_fit_rep$summary.nu$mode
    )
  )
  print(result_bru_rep_df)

##   parameter true      mean      mode
## 1   std.dev  1.3 1.2808406 1.2801461
## 2     range  0.2 0.1912403 0.1937536
## 3        nu  0.9 0.8970910 0.8602582

We can also plot the posterior marginal densities with the help of the gg_df() function:

  posterior_df_bru_fit_rep <- gg_df(result_bru_fit_rep)

  ggplot(posterior_df_bru_fit_rep) + geom_line(aes(x = x, y = y)) + 
  facet_wrap(~parameter, scales = "free") + labs(y = "Density")